; ===========================================================================;
; Data: surface pressure
; Grid: Lat/long grid
; Data Format: netCDF
; drs 03/10/2003
; ===========================================================================;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load "./functions_contrib.ncl"
;load "./functions_surfaces.ncl"
;load "./funcs_surf_variance.ncl"
; ===========================================================================;

begin
seas=(/"DJF","JJA","MAM","SON","ANN"/)
nseas=(/"DJF","JJA","MAM","SON","ANN"/)
reg=(/"india_eafrica","africa","global"/)

do ir=1,1	;region
do i=4,4	;season
  outputfile = "ts_seas_" + seas(i) +"_1970-1989_diff_" + reg(ir)

;  plot = new(6,graphic)                          ; create a plot array
;  dum = new(6,graphic)                          ; create a plot array
;  dum1 = new(6,graphic)                          ; create a plot array
 ; dum2 = new(6,graphic)                          ; create a plot array
;  dum3 = new(6,graphic)                          ; create a plot array
;  dum4 = new(6,graphic)                          ; create a plot array

 inputfile1 = "/home/divanova/SahelProject/amwg/climo/e.e122.E_1850_CAM5.f09_g16.t4Ind_stg_co2b_cam5.20-49y/e.e122.E_1850_CAM5.f09_g16.t4Ind_stg_co2b_cam5_" + seas(i) +"_climo.nc"
 inputfile2 = "/home/divanova/SahelProject/amwg/climo/f.e122.F20C5TR.f09_g16.t4anom.70-89y/f.e122.F20C5TR.f09_g16.t4anom_" + seas(i) +"_climo.nc"
 inputfile3 = "/home/divanova/SahelProject/amwg/climo/e.e122.E_1850_CAM5.f09_g16.sahel_controlCO2III_cam5.20-49y/e.e122.E_1850_CAM5.f09_g16.sahel_controlCO2III_cam5_" + seas(i) +"_climo.nc"
 inputfile4 = "/home/divanova/SahelProject/amwg/climo/f.e122.F20C5TR.f09_g16.1970-1989y/f.e122.F20C5TR.f09_g16_" + seas(i) + "_climo.nc"

;        shape files


 dir="/home/divanova/SahelProject/topo/"
 filename="sahel_zone/sahelzone.shp"
 filename1="India/India_Boundary_Igismap/India_Boundary.shp"
 filename2="TANZANIA/gadm36_TZA_0.shp"
 filename3="KENYA/Igismap_Kenya/Kenya_Boundary.shp"
 filename4="UGANDA/gadm36_UGA_0.shp"

  SigLvl=0.05

  f1 = addfile(inputfile1,"r")
  f2 = addfile(inputfile2,"r")
  f3 = addfile(inputfile3,"r")
  f4 = addfile(inputfile4,"r")

  ts1=f1->TS(0,:,:)
  m1=f1->OCNFRAC(0,:,:)
  u1=ts1
  u1=mask(ts1,m1,1)
  u2=f2->TS(0,:,:)
  ts3=f3->TS(0,:,:)
  m3=f3->OCNFRAC(0,:,:)
  u3=ts3
  u3=mask(ts3,m3,1)
  u4=f4->TS(0,:,:)
  u1=u1-273.15
  u2=u2-273.15
  u3=u3-273.15
  u4=u4-273.15

  diff_u1=u1
  diff_u1=u1-u3
  diff_u2=u2
  diff_u2=u2-u4

  wks  = gsn_open_wks("eps",outputfile)        ; open a ps file
  print(outputfile)
; wks  = gsn_open_wks("X11","")               ; open an X11 window file
;
;  gsn_define_colormap(wks,"BlWhRe")      ; choose coloramp
;  gsn_define_colormap(wks,"gui_default")      ; choose coloramp
;  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")      ; choose coloramp
;  gsn_define_colormap(wks,"cmocean_ice")      ; choose coloramp
;  gsn_define_colormap(wks,"WhiteBlue")      ; choose coloramp
;  gsn_reverse_colormap(wks) 
;  setvalues wks            
;    "wkForegroundColor" : (/1.,1.,1./) 
;    "wkBackgroundColor" : (/1.,1.,1./)  
;    "wkBackgroundColor" : (/0.,0.,0./)  
;  end setvalues
;  gsn_define_colormap(wks,"temp_diff_18lev")      ; choose coloramp
  gsn_define_colormap(wks,"cmp_b2r")      ; choose coloramp
;  gsn_define_colormap(wks,"ncl_default")      ; choose coloramp

; =================================================;
; create vector plot
; =================================================;
  res                   = True                ; plot mods desired
  res@gsnDraw              = False               ; do not draw picture
  res@gsnFrame             = False               ; do not advance frame
  res@trGridType          = "TriangularMesh"
  res@gsnAddCyclic = True
  res@mpPerimOn   = True        ; Turn on map perimeter.

  res@mpFillOn              = False              ; turn off map fill
  res@mpOutlineOn           = True               ; turn on map outlines
  res@mpCenterLonF      = 0                ; center longitude

  if (reg(ir).eq."global") then
  res@mpMinLonF         = -180                  ; minimum lat to plot
  res@mpMaxLonF         = 180                  ; minimum lat to plot
  res@mpMinLatF         = -90                  ; minimum lat to plot
  res@mpMaxLatF         = 90                  ; minimum lat to plot
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  -25             ; set min contour level
  res@cnMaxLevelValF       =  40             ; set max contour level
  res@cnLevelSpacingF   = 5                   ; interval spacing
  end if

  if (reg(ir).eq."india_eafrica") then
  res@mpMinLonF         = -30                  ; minimum lat to plot
  res@mpMaxLonF         = 120                  ; minimum lat to plot
  res@mpMinLatF         = -60                  ; minimum lat to plot
  res@mpMaxLatF         = 60                  ; minimum lat to plot
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  -25             ; set min contour level
  res@cnMaxLevelValF       =  40             ; set max contour level
  res@cnLevelSpacingF   = 5                   ; interval spacing
  end if

; Africa
  if (reg(ir).eq."africa") then
  res@mpMinLonF         = -20                  ; minimum lat to plot
  res@mpMaxLonF         = 60                  ; minimum lat to plot
  res@mpMinLatF         = -30                  ; minimum lat to plot
  res@mpMaxLatF         = 30                  ; minimum lat to plot
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  10             ; set min contour level
  res@cnMaxLevelValF       =  40             ; set max contour level
  res@cnLevelSpacingF   = 5                   ; interval spacing
  end if


  res@cnFillOn          = True                  ; color fill
; when RasterFill is on only Solid FillPatern is  possible
;  res@cnFillMode = "RasterFill"
;  res@cnRasterSmoothingOn = True

  res@cnLinesOn         = False                 ; no contour lines
  res@cnLineLabelsOn         = False                 ; no contour lines

;  res@cnMissingValFillColor="white" ; filled values drawn in gray
;;  res@cnMissingValFillColor="gray70" ; filled values drawn in gray
;  res@cnMissingValFillPattern=17 ; filled valyes are dotted
;  res@cnFillBackgroundColor=-1 ; contour background color is transparent


  res@gsnSpreadColors     = True              ; use full colormap
  res@gsnSpreadColorEnd   = -3                ; don't use last colors (grey)
  res@lbLabelBarOn        = True               ; Turns off the individual label bars

  res@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
  res@gsnMajorLonSpacing = 30              ; change maj lon tm spacing
  res@tmXBMinorOn        = False           ; no lon minor tickmarks
  res@tmYLMinorOn        = False           ; no lat minor tickmarks
  res@tmYLOn       = False           ; no lat minor tickmarks
  res@tmXTOn       = False           ; no lat minor tickmarks


  plres = True
  plres@gsLineColor = "magenta"
  plres@gsLineThicknessF =5

  plres1 = True
  plres1@gsLineColor = "green"
  plres1@gsLineThicknessF =2


  res1 = True
;  res1@gsnDraw              = False               ; do not draw picture
;  res1@gsnFrame             = False               ; do not advance frame
  res1@mpFillOn              = True              ; turn off map fill
  res1@mpLandFillColor         = "grey"          ; make land white

  res1@mpOutlineOn           = True               ; turn on map outlines
  res1@mpCenterLonF      = 0                ; center longitude
  res1@mpPerimOn   = True        ; Turn on map perimeter.


  res1@gsnMajorLatSpacing = 30              ; change maj lat tm spacing
  res1@gsnMajorLonSpacing = 30              ; change maj lon tm spacing
  res1@tmXBMinorOn        = False           ; no lon minor tickmarks
  res1@tmYLMinorOn        = False           ; no lat minor tickmarks
  res1@tmYLOn       = False           ; no lat minor tickmarks
  res1@tmXTOn       = False           ; no lat minor tickmarks

  if (reg(ir).eq."global") then
  res1@mpMinLonF         = -180                  ; minimum lat to plot
  res1@mpMaxLonF         = 180                  ; minimum lat to plot
  res1@mpMinLatF         = -90                  ; minimum lat to plot
  res1@mpMaxLatF         = 90                  ; minimum lat to plot
  end if

  if (reg(ir).eq."india_eafrica") then
  res1@mpMinLonF         = -30                  ; minimum lat to plot
  res1@mpMaxLonF         = 120                  ; minimum lat to plot
  res1@mpMinLatF         = -60                  ; minimum lat to plot
  res1@mpMaxLatF         = 60                  ; minimum lat to plot
  end if

; Africa
  if (reg(ir).eq."africa") then
  res1@mpMinLonF         = -20                  ; minimum lat to plot
  res1@mpMaxLonF         = 60                  ; minimum lat to plot
  res1@mpMinLatF         = -30                  ; minimum lat to plot
  res1@mpMaxLatF         = 30                  ; minimum lat to plot
  end if

  res1@gsnAddCyclic = True
  res1@cnFillPalette     = "BlueDarkRed18"
  res1@cnFillOn          = True                  ; color fill
  res1@cnLinesOn         = False                 ; no contour lines
  res1@cnLineLabelsOn         = False                 ; no contour lines
  res1@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res1@cnMinLevelValF       =  -3.             ; set min contour level
  res1@cnMaxLevelValF       =  3.             ; set max contour level
  res1@cnLevelSpacingF   = .2                   ; interval spacing
  res1@mpLandFillColor         = "grey"          ; make land white


  res1@gsnRightString          = nseas(i) + ", 20-49y" 
  res1@gsnLeftString          = "io strongly cooled-2xCO2"   ; title  
  plot = gsn_csm_contour_map(wks,diff_u1,res1)
  dum  = gsn_add_shapefile_polylines(wks,plot,dir+filename,plres)
;  dum1  = gsn_add_shapefile_polylines(wks,plot,dir+filename1,plres1)
;  dum2  = gsn_add_shapefile_polylines(wks,plot,dir+filename2,plres1)
;  dum3(4)  = gsn_add_shapefile_polylines(wks,plot(4),dir+filename3,plres1)
;  dum4(4)  = gsn_add_shapefile_polylines(wks,plot(4),dir+filename4,plres1)




end do
end do
end

